home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
basic
/
chaosexe.zip
/
XPOINCAR.TRU
< prev
next >
Wrap
Text File
|
1980-01-01
|
3KB
|
113 lines
!PROGRAM TITLE "XPOINCAR(E)"
CLEAR
PRINT" ***PENDULUM - POINCARE SECTION***"
PRINT
PRINT"THIS PROGRAM DISPLAYS THE POINCARE SECTION OF THE PENDULUM"
PRINT"AND CAN SAVE THE DATA TO A FILE."
LIBRARY "SGLIB.TRC"
!
DECLARE DEF accel
DIM A(1),B(1)
INPUT prompt"Input driving force strength: ":g
INPUT prompt"Input damping (If no damping then input 9999999):":q
INPUT prompt"Input initial angle, angular velocity: ": xint,vint
INPUT Prompt"Input min. and max. time:":tmin,tmax
INPUT prompt"Input phase angle/(2*pi): ":phi
INPUT PROMPT" SAVE DATA TO A FILE? YES(1), NO(2):":SAVEFILE
IF SAVEFILE=1 THEN
INPUT PROMPT"FILE NAME FORMAT EX. 14954020 :":FILENAME
INPUT PROMPT"DRIVE FOR FILE DISK A,B,C,ETC.:":DISK$
LET NAME$=DISK$&":"&STR$(FILENAME)
END IF
!
CALL PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)
CALL SETXSCALE(XMIN,XMAX)
CALL SETYSCALE(YMIN,YMAX)
CALL SETTEXT("PENDULUM POINCARE SECTION","THETA","OMEGA")
CALL RESERVELEGEND
!
DATA O,O
CALL DATAGRAPH(A,B,1,O,"RED")
LET t=0
LET x=xint
LET v=vint
CALL GOTOCANVAS
!
!CALCULATION AND GRAPHING BLOCK
LET phi=phi*2*pi
IF SAVEFILE=1 THEN
OPEN #1:NAME NAME$, ORGANIZATION RECORD, CREATE NEWOLD
ASK #1:FILESIZE LENGTH
IF LENGTH=0 THEN SET#1:RECSIZE 10
SET #1: POINTER END
END IF
FOR i=1 to 1000000
CALL rk4(x,v,tstep,xnew,vnew,t,w,g,q) ! Take a step of size tstep
LET tshalf=tstep/2
CALL rk4(x,v,tshalf,xnh,vnh,t,w,g,q) !Take two half steps
CALL rk4(xnh,vnh,tshalf,xn,vn,t+tshalf,w,g,q)
LET d1=abs(xn-xnew)
LET d2=abs(vn-vnew)
LET delta=max(d1,d2)
IF delta<eps then
IF t>tmin then
LET tnew=t+tstep
LET w1=mod(phi-w*t,2*pi) !Check for Poincare section
LET w2=mod(w*tnew-phi,2*pi)
IF w1<w*tstep then
IF w2<w*tstep then
LET ts=w1/w
CALL rk4(x,v,ts,xp,vp,t,w,g,q) !CALCULATES POINT AT SECTION
IF abs(xp)>pi then LET xp=xp-2*pi*abs(xp)/xp
CALL GRAPHPOINT(XP,VP,1)
IF SAVEFILE=1 THEN WRITE #1:XP,VP
END IF
END IF
END IF
LET x=xnew
LET v=vnew
LET t=t+tstep !Expand step size
LET tstep=tstep*.95*(eps/delta)^.25
IF abs(x)>pi then !bring theta back into range
LET x=x-2*pi*abs(x)/x
END IF
ELSE !else reduce step size
LET tstep=tstep*.95*(eps/delta)^.2
END IF
IF t>tmax then LET i=1000001
NEXT i
LET G$=STR$(G)
LET Q$=STR$(Q)
CALL ADDLEGEND("G="&STR$(G)&" Q="&STR$(Q)&" PHI="&STR$(PHI),0,1,"WHITE")
CALL DRAWLEGEND
END
!
!
SUB rk4(x,v,tstep,xnew,vnew,t,w,g,q)
DECLARE DEF accel
LET xk1=tstep*v
LET vk1=tstep*accel(x,v,t,w,g,q)
LET xk2=tstep*(v+vk1/2)
LET vk2=tstep*accel(x+xk1/2,v+vk1/2,t+tstep/2,w,g,q)
LET xk3=tstep*(v+vk2/2)
LET vk3=tstep*accel(x+xk2/2,v+vk2/2,t+tstep/2,w,g,q)
LET xk4=tstep*(v+vk3)
LET vk4=tstep*accel(x+xk3,v+vk3,t+tstep,w,g,q)
LET vnew=v+(vk1+2*vk2+2*vk3+vk4)/6
LET xnew=x+(xk1+2*xk2+2*xk3+xk4)/6
END SUB
DEF accel(x,v,t,w,g,q)
LET accel=-sin(x)-(1/q)*v+g*cos(w*t)
END def
!
SUB PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)
LET W=0.66666666
LET EPS=1.0E-6
LET TSTEP=0.5
LET XMIN=-3
LET XMAX=3
LET YMIN=-3
LET YMAX=3
END SUB